//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
///////////////////////////////////////////////////////////////////////////////
//
// MODULE:        ExSPSPosDistribution.cc
//
// Version:      1.0
// Date:         5/02/04
// Author:       Fan Lei
// Organisation: QinetiQ ltd.
// Customer:     ESA/ESTEC
//
///////////////////////////////////////////////////////////////////////////////
//
// CHANGE HISTORY
// --------------
//
//
// Version 1.0, 05/02/2004, Fan Lei, Created.
//    Based on the G4GeneralParticleSource class in Geant4 v6.0
//
// 2016, Linyan, Modified to include CopyFromSPSP.
//
///////////////////////////////////////////////////////////////////////////////
//

#include "ExSPSPosDistribution.hh"

#include "G4PhysicalConstants.hh"
#include "Randomize.hh"
#include "G4TransportationManager.hh"
#include "G4VPhysicalVolume.hh"
#include "G4PhysicalVolumeStore.hh"

ExSPSPosDistribution::ExSPSPosDistribution()
    : posRndm(0)
{

  // Initialise all variables
  // Position distribution Variables

  SourcePosType = "Point";
  Shape = "NULL";
  halfx = 0.;
  halfy = 0.;
  halfz = 0.;
  Radius = 0.;
  Radius0 = 0.;
  SR = 0.;
  SX = 0.;
  SY = 0.;
  ParAlpha = 0.;
  ParTheta = 0.;
  ParPhi = 0.;
  CentreCoords = G4ThreeVector(0., 0., 0.);
  Rotx = CLHEP::HepXHat;
  Roty = CLHEP::HepYHat;
  Rotz = CLHEP::HepZHat;
  Confine = false; // If true confines source distribution to VolName
  VolName = "NULL";
  SideRefVec1 = CLHEP::HepXHat; // x-axis
  SideRefVec2 = CLHEP::HepYHat; // y-axis
  SideRefVec3 = CLHEP::HepZHat; // z-axis
  verbosityLevel = 0;
  gNavigator = G4TransportationManager::GetTransportationManager()
                   ->GetNavigatorForTracking();
}

ExSPSPosDistribution::~ExSPSPosDistribution()
{
}

void ExSPSPosDistribution::SetPosDisType(G4String PosType)
{
  SourcePosType = PosType;
}

void ExSPSPosDistribution::SetPosDisShape(G4String shapeType)
{
  Shape = shapeType;
}

void ExSPSPosDistribution::SetCentreCoords(G4ThreeVector coordsOfCentre)
{
  CentreCoords = coordsOfCentre;
}

void ExSPSPosDistribution::SetPosRot1(G4ThreeVector posrot1)
{
  // This should be x'
  Rotx = posrot1;
  if (verbosityLevel == 2)
  {
    G4cout << "Vector x' " << Rotx << G4endl;
  }
  GenerateRotationMatrices();
}

void ExSPSPosDistribution::SetPosRot2(G4ThreeVector posrot2)
{
  // This is a vector in the plane x'y' but need not
  // be y'
  Roty = posrot2;
  if (verbosityLevel == 2)
  {
    G4cout << "The vector in the x'-y' plane " << Roty << G4endl;
  }
  GenerateRotationMatrices();
}

void ExSPSPosDistribution::SetHalfX(G4double xhalf)
{
  halfx = xhalf;
}

void ExSPSPosDistribution::SetHalfY(G4double yhalf)
{
  halfy = yhalf;
}

void ExSPSPosDistribution::SetHalfZ(G4double zhalf)
{
  halfz = zhalf;
}

void ExSPSPosDistribution::SetRadius(G4double rds)
{
  Radius = rds;
}

void ExSPSPosDistribution::SetRadius0(G4double rds)
{
  Radius0 = rds;
}

void ExSPSPosDistribution::SetBeamSigmaInR(G4double r)
{
  SX = SY = r;
  SR = r;
}

void ExSPSPosDistribution::SetBeamSigmaInX(G4double r)
{
  SX = r;
}

void ExSPSPosDistribution::SetBeamSigmaInY(G4double r)
{
  SY = r;
}

void ExSPSPosDistribution::SetParAlpha(G4double paralp)
{
  ParAlpha = paralp;
}

void ExSPSPosDistribution::SetParTheta(G4double parthe)
{
  ParTheta = parthe;
}

void ExSPSPosDistribution::SetParPhi(G4double parphi)
{
  ParPhi = parphi;
}

void ExSPSPosDistribution::GenerateRotationMatrices()
{
  // This takes in 2 vectors, x' and one in the plane x'-y',
  // and from these takes a cross product to calculate z'.
  // Then a cross product is taken between x' and z' to give
  // y'.
  Rotx = Rotx.unit();      // x'
  Roty = Roty.unit();      // vector in x'y' plane
  Rotz = Rotx.cross(Roty); // z'
  Rotz = Rotz.unit();
  Roty = Rotz.cross(Rotx); // y'
  Roty = Roty.unit();
  if (verbosityLevel == 2)
  {
    G4cout << "The new axes, x', y', z' " << Rotx << " " << Roty << " " << Rotz << G4endl;
  }
}

void ExSPSPosDistribution::ConfineSourceToVolume(G4String Vname)
{
  VolName = Vname;
  if (verbosityLevel == 2)
    G4cout << VolName << G4endl;
  G4VPhysicalVolume *tempPV = NULL;
  G4PhysicalVolumeStore *PVStore = 0;
  G4String theRequiredVolumeName = VolName;
  PVStore = G4PhysicalVolumeStore::GetInstance();
  G4int i = 0;
  G4bool found = false;
  if (verbosityLevel == 2)
    G4cout << PVStore->size() << G4endl;
  while (!found && i < G4int(PVStore->size()))
  {
    tempPV = (*PVStore)[i];
    found = tempPV->GetName() == theRequiredVolumeName;
    if (verbosityLevel == 2)
      G4cout << i << " "
             << " " << tempPV->GetName() << " " << theRequiredVolumeName << " " << found << G4endl;
    if (!found)
    {
      i++;
    }
  }
  // found = true then the volume exists else it doesnt.
  if (found == true)
  {
    if (verbosityLevel >= 1)
      G4cout << "Volume " << VolName << " exists" << G4endl;
    Confine = true;
  }
  else if (theRequiredVolumeName != "NULL")
  {
    G4cout << " **** Error: Volume does not exist **** " << G4endl;
    G4cout << theRequiredVolumeName << G4endl;
    G4cout << " Ignoring confine condition" << G4endl;
    Confine = false;
    VolName = "NULL";
  }
}

void ExSPSPosDistribution::GeneratePointSource()
{
  // Generates Points given the point source.
  if (SourcePosType == "Point")
    particle_position = CentreCoords;
  else if (verbosityLevel >= 1)
    G4cout << "Error SourcePosType is not set to Point" << G4endl;
}

void ExSPSPosDistribution::GeneratePointsInBeam()
{
  G4double x, y, z;

  G4ThreeVector RandPos;
  G4double tempx, tempy, tempz;
  z = 0.;

  // Private Method to create points in a plane
  if (Shape == "Circle")
  {
    x = Radius + 100.;
    y = Radius + 100.;
    while (std::sqrt((x * x) + (y * y)) > Radius)
    {
      x = posRndm->GenRandX();
      y = posRndm->GenRandY();

      x = (x * 2. * Radius) - Radius;
      y = (y * 2. * Radius) - Radius;
    }
    x += G4RandGauss::shoot(0.0, SX);
    y += G4RandGauss::shoot(0.0, SY);
  }
  else
  {
    // all other cases default to Rectangle case
    x = posRndm->GenRandX();
    y = posRndm->GenRandY();
    x = (x * 2. * halfx) - halfx;
    y = (y * 2. * halfy) - halfy;
    x += G4RandGauss::shoot(0.0, SX);
    y += G4RandGauss::shoot(0.0, SY);
  }
  // Apply Rotation Matrix
  // x * Rotx, y * Roty and z * Rotz
  if (verbosityLevel >= 2)
  {
    G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
  }
  tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
  tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
  tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());

  RandPos.setX(tempx);
  RandPos.setY(tempy);
  RandPos.setZ(tempz);

  // Translate
  particle_position = CentreCoords + RandPos;
  if (verbosityLevel >= 1)
  {
    if (verbosityLevel >= 2)
    {
      G4cout << "Rotated Position " << RandPos << G4endl;
    }
    G4cout << "Rotated and Translated position " << particle_position << G4endl;
  }
}

void ExSPSPosDistribution::GeneratePointsInPlane()
{
  G4double x, y, z;
  G4double expression;
  G4ThreeVector RandPos;
  G4double tempx, tempy, tempz;
  x = y = z = 0.;

  if (SourcePosType != "Plane" && verbosityLevel >= 1)
    G4cout << "Error: SourcePosType is not Plane" << G4endl;

  // Private Method to create points in a plane
  if (Shape == "Circle")
  {
    x = Radius + 100.;
    y = Radius + 100.;
    while (std::sqrt((x * x) + (y * y)) > Radius)
    {
      x = posRndm->GenRandX();
      y = posRndm->GenRandY();

      x = (x * 2. * Radius) - Radius;
      y = (y * 2. * Radius) - Radius;
    }
  }
  else if (Shape == "Annulus")
  {
    x = Radius + 100.;
    y = Radius + 100.;
    while (std::sqrt((x * x) + (y * y)) > Radius || std::sqrt((x * x) + (y * y)) < Radius0)
    {
      x = posRndm->GenRandX();
      y = posRndm->GenRandY();

      x = (x * 2. * Radius) - Radius;
      y = (y * 2. * Radius) - Radius;
    }
  }
  else if (Shape == "Ellipse")
  {
    expression = 20.;
    while (expression > 1.)
    {
      x = posRndm->GenRandX();
      y = posRndm->GenRandY();

      x = (x * 2. * halfx) - halfx;
      y = (y * 2. * halfy) - halfy;

      expression = ((x * x) / (halfx * halfx)) + ((y * y) / (halfy * halfy));
    }
  }
  else if (Shape == "Square")
  {
    x = posRndm->GenRandX();
    y = posRndm->GenRandY();
    x = (x * 2. * halfx) - halfx;
    y = (y * 2. * halfy) - halfy;
  }
  else if (Shape == "Rectangle")
  {
    x = posRndm->GenRandX();
    y = posRndm->GenRandY();
    x = (x * 2. * halfx) - halfx;
    y = (y * 2. * halfy) - halfy;
  }
  else
    G4cout << "Shape not one of the plane types" << G4endl;

  // Apply Rotation Matrix
  // x * Rotx, y * Roty and z * Rotz
  if (verbosityLevel == 2)
  {
    G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
  }
  tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
  tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
  tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());

  RandPos.setX(tempx);
  RandPos.setY(tempy);
  RandPos.setZ(tempz);

  // Translate
  particle_position = CentreCoords + RandPos;
  if (verbosityLevel >= 1)
  {
    if (verbosityLevel == 2)
    {
      G4cout << "Rotated Position " << RandPos << G4endl;
    }
    G4cout << "Rotated and Translated position " << particle_position << G4endl;
  }

  // For Cosine-Law make SideRefVecs = to Rotation matrix vectors
  SideRefVec1 = Rotx;
  SideRefVec2 = Roty;
  SideRefVec3 = Rotz;
  // If rotation matrix z' point to origin then invert the matrix
  // So that SideRefVecs point away.
  if ((CentreCoords.x() > 0. && Rotz.x() < 0.) || (CentreCoords.x() < 0. && Rotz.x() > 0.) || (CentreCoords.y() > 0. && Rotz.y() < 0.) || (CentreCoords.y() < 0. && Rotz.y() > 0.) || (CentreCoords.z() > 0. && Rotz.z() < 0.) || (CentreCoords.z() < 0. && Rotz.z() > 0.))
  {
    // Invert y and z.
    SideRefVec2 = -SideRefVec2;
    SideRefVec3 = -SideRefVec3;
  }
  if (verbosityLevel == 2)
  {
    G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
  }
}

void ExSPSPosDistribution::GeneratePointsOnSurface()
{
  // Private method to create points on a surface
  G4double theta, phi;
  G4double x, y, z;
  x = y = z = 0.;
  G4ThreeVector RandPos;
  //  G4double tempx, tempy, tempz;

  if (SourcePosType != "Surface" && verbosityLevel >= 1)
    G4cout << "Error SourcePosType not Surface" << G4endl;

  if (Shape == "Sphere")
  {
    // G4double tantheta;
    theta = posRndm->GenRandPosTheta();
    phi = posRndm->GenRandPosPhi();
    theta = std::acos(1. - 2. * theta); // theta isotropic
    phi = phi * 2. * pi;
    // tantheta = std::tan(theta);

    x = Radius * std::sin(theta) * std::cos(phi);
    y = Radius * std::sin(theta) * std::sin(phi);
    z = Radius * std::cos(theta);

    RandPos.setX(x);
    RandPos.setY(y);
    RandPos.setZ(z);

    // Cosine-law (not a good idea to use this here)
    G4ThreeVector zdash(x, y, z);
    zdash = zdash.unit();
    G4ThreeVector xdash = Rotz.cross(zdash);
    G4ThreeVector ydash = xdash.cross(zdash);
    SideRefVec1 = xdash.unit();
    SideRefVec2 = ydash.unit();
    SideRefVec3 = zdash.unit();
  }
  else if (Shape == "Ellipsoid")
  {
    G4double minphi, maxphi, middlephi;
    G4double answer, constant;

    constant = pi / (halfx * halfx) + pi / (halfy * halfy) +
               twopi / (halfz * halfz);

    // simplified approach
    theta = posRndm->GenRandPosTheta();
    phi = posRndm->GenRandPosPhi();

    theta = std::acos(1. - 2. * theta);
    minphi = 0.;
    maxphi = twopi;
    while (maxphi - minphi > 0.)
    {
      middlephi = (maxphi + minphi) / 2.;
      answer = (1. / (halfx * halfx)) * (middlephi / 2. + std::sin(2 * middlephi) / 4.) + (1. / (halfy * halfy)) * (middlephi / 2. - std::sin(2 * middlephi) / 4.) + middlephi / (halfz * halfz);
      answer = answer / constant;
      if (answer > phi)
        maxphi = middlephi;
      if (answer < phi)
        minphi = middlephi;
      if (std::fabs(answer - phi) <= 0.00001)
      {
        minphi = maxphi + 1;
        phi = middlephi;
      }
    }

    x = std::sin(theta) * std::cos(phi);
    y = std::sin(theta) * std::sin(phi);
    z = std::cos(theta);
    // x,y and z form a unit vector. Put this onto the ellipse.
    G4double lhs;
    // solve for x
    G4double numYinX = y / x;
    G4double numZinX = z / x;
    G4double tempxvar;
    tempxvar = 1. / (halfx * halfx) + (numYinX * numYinX) / (halfy * halfy) + (numZinX * numZinX) / (halfz * halfz);

    tempxvar = 1. / tempxvar;
    G4double coordx = std::sqrt(tempxvar);

    // solve for y
    G4double numXinY = x / y;
    G4double numZinY = z / y;
    G4double tempyvar;
    tempyvar = (numXinY * numXinY) / (halfx * halfx) + 1. / (halfy * halfy) + (numZinY * numZinY) / (halfz * halfz);
    tempyvar = 1. / tempyvar;
    G4double coordy = std::sqrt(tempyvar);

    // solve for z
    G4double numXinZ = x / z;
    G4double numYinZ = y / z;
    G4double tempzvar;
    tempzvar = (numXinZ * numXinZ) / (halfx * halfx) + (numYinZ * numYinZ) / (halfy * halfy) + 1. / (halfz * halfz);
    tempzvar = 1. / tempzvar;
    G4double coordz = std::sqrt(tempzvar);

    lhs = std::sqrt((coordx * coordx) / (halfx * halfx) +
                    (coordy * coordy) / (halfy * halfy) +
                    (coordz * coordz) / (halfz * halfz));

    if (std::fabs(lhs - 1.) > 0.001 && verbosityLevel >= 1)
      G4cout << "Error: theta, phi not really on ellipsoid" << G4endl;

    // coordx, coordy and coordz are all positive
    G4double TestRandVar = G4UniformRand();
    if (TestRandVar > 0.5)
    {
      coordx = -coordx;
    }
    TestRandVar = G4UniformRand();
    if (TestRandVar > 0.5)
    {
      coordy = -coordy;
    }
    TestRandVar = G4UniformRand();
    if (TestRandVar > 0.5)
    {
      coordz = -coordz;
    }

    RandPos.setX(coordx);
    RandPos.setY(coordy);
    RandPos.setZ(coordz);

    // Cosine-law (not a good idea to use this here)
    G4ThreeVector zdash(coordx, coordy, coordz);
    zdash = zdash.unit();
    G4ThreeVector xdash = Rotz.cross(zdash);
    G4ThreeVector ydash = xdash.cross(zdash);
    SideRefVec1 = xdash.unit();
    SideRefVec2 = ydash.unit();
    SideRefVec3 = zdash.unit();
  }
  else if (Shape == "Cylinder")
  {
    G4double AreaTop, AreaBot, AreaLat;
    G4double AreaTotal, prob1, prob2, prob3;
    G4double testrand;

    // User giver Radius and z-half length
    // Calculate surface areas, maybe move this to
    // a different routine.

    AreaTop = pi * Radius * Radius;
    AreaBot = AreaTop;
    AreaLat = 2. * pi * Radius * 2. * halfz;
    AreaTotal = AreaTop + AreaBot + AreaLat;

    prob1 = AreaTop / AreaTotal;
    prob2 = AreaBot / AreaTotal;
    prob3 = 1.00 - prob1 - prob2;
    if (std::fabs(prob3 - (AreaLat / AreaTotal)) >= 0.001)
    {
      if (verbosityLevel >= 1)
        G4cout << AreaLat / AreaTotal << " " << prob3 << G4endl;
      G4cout << "Error in prob3" << G4endl;
    }

    // Decide surface to calculate point on.

    testrand = G4UniformRand();
    if (testrand <= prob1)
    {
      // Point on Top surface
      z = halfz;
      x = Radius + 100.;
      y = Radius + 100.;
      while (((x * x) + (y * y)) > (Radius * Radius))
      {
        x = posRndm->GenRandX();
        y = posRndm->GenRandY();

        x = x * 2. * Radius;
        y = y * 2. * Radius;
        x = x - Radius;
        y = y - Radius;
      }
      // Cosine law
      SideRefVec1 = Rotx;
      SideRefVec2 = Roty;
      SideRefVec3 = Rotz;
    }
    else if ((testrand > prob1) && (testrand <= (prob1 + prob2)))
    {
      // Point on Bottom surface
      z = -halfz;
      x = Radius + 100.;
      y = Radius + 100.;
      while (((x * x) + (y * y)) > (Radius * Radius))
      {
        x = posRndm->GenRandX();
        y = posRndm->GenRandY();

        x = x * 2. * Radius;
        y = y * 2. * Radius;
        x = x - Radius;
        y = y - Radius;
      }
      // Cosine law
      SideRefVec1 = Rotx;
      SideRefVec2 = -Roty;
      SideRefVec3 = -Rotz;
    }
    else if (testrand > (prob1 + prob2))
    {
      G4double rand;
      // Point on Lateral Surface

      rand = posRndm->GenRandPosPhi();
      rand = rand * 2. * pi;

      x = Radius * std::cos(rand);
      y = Radius * std::sin(rand);

      z = posRndm->GenRandZ();

      z = z * 2. * halfz;
      z = z - halfz;

      // Cosine law
      G4ThreeVector zdash(x, y, 0.);
      zdash = zdash.unit();
      G4ThreeVector xdash = Rotz.cross(zdash);
      G4ThreeVector ydash = xdash.cross(zdash);
      SideRefVec1 = xdash.unit();
      SideRefVec2 = ydash.unit();
      SideRefVec3 = zdash.unit();
    }
    else
      G4cout << "Error: testrand " << testrand << G4endl;

    RandPos.setX(x);
    RandPos.setY(y);
    RandPos.setZ(z);
  }
  else if (Shape == "Para")
  {
    G4double testrand;
    // Right Parallelepiped.
    //  User gives x,y,z half lengths and ParAlpha
    //  ParTheta and ParPhi
    //  +x = <1, -x >1 & <2, +y >2 & <3, -y >3 &<4
    //  +z >4 & < 5, -z >5 &<6.
    testrand = G4UniformRand();
    G4double AreaX = halfy * halfz * 4.;
    G4double AreaY = halfx * halfz * 4.;
    G4double AreaZ = halfx * halfy * 4.;
    G4double AreaTotal = 2 * (AreaX + AreaY + AreaZ);
    G4double Probs[6];
    Probs[0] = AreaX / AreaTotal;
    Probs[1] = Probs[0] + AreaX / AreaTotal;
    Probs[2] = Probs[1] + AreaY / AreaTotal;
    Probs[3] = Probs[2] + AreaY / AreaTotal;
    Probs[4] = Probs[3] + AreaZ / AreaTotal;
    Probs[5] = Probs[4] + AreaZ / AreaTotal;

    x = posRndm->GenRandX();
    y = posRndm->GenRandY();
    z = posRndm->GenRandZ();

    x = x * halfx * 2.;
    x = x - halfx;
    y = y * halfy * 2.;
    y = y - halfy;
    z = z * halfz * 2.;
    z = z - halfz;
    // Pick a side first
    if (testrand < Probs[0])
    {
      // side is +x
      x = halfx + z * std::tan(ParTheta) * std::cos(ParPhi) + y * std::tan(ParAlpha);
      y = y + z * std::tan(ParTheta) * std::sin(ParPhi);
      // z = z;
      // Cosine-law
      G4ThreeVector xdash(halfz * std::tan(ParTheta) * std::cos(ParPhi),
                          halfz * std::tan(ParTheta) * std::sin(ParPhi),
                          halfz / std::cos(ParPhi));
      G4ThreeVector ydash(halfy * std::tan(ParAlpha), -halfy, 0.0);
      xdash = xdash.unit();
      ydash = ydash.unit();
      G4ThreeVector zdash = xdash.cross(ydash);
      SideRefVec1 = xdash.unit();
      SideRefVec2 = ydash.unit();
      SideRefVec3 = zdash.unit();
    }
    else if (testrand >= Probs[0] && testrand < Probs[1])
    {
      // side is -x
      x = -halfx + z * std::tan(ParTheta) * std::cos(ParPhi) + y * std::tan(ParAlpha);
      y = y + z * std::tan(ParTheta) * std::sin(ParPhi);
      // z = z;
      // Cosine-law
      G4ThreeVector xdash(halfz * std::tan(ParTheta) * std::cos(ParPhi),
                          halfz * std::tan(ParTheta) * std::sin(ParPhi),
                          halfz / std::cos(ParPhi));
      G4ThreeVector ydash(halfy * std::tan(ParAlpha), halfy, 0.0);
      xdash = xdash.unit();
      ydash = ydash.unit();
      G4ThreeVector zdash = xdash.cross(ydash);
      SideRefVec1 = xdash.unit();
      SideRefVec2 = ydash.unit();
      SideRefVec3 = zdash.unit();
    }
    else if (testrand >= Probs[1] && testrand < Probs[2])
    {
      // side is +y
      x = x + z * std::tan(ParTheta) * std::cos(ParPhi) + halfy * std::tan(ParAlpha);
      y = halfy + z * std::tan(ParTheta) * std::sin(ParPhi);
      // z = z;
      // Cosine-law
      G4ThreeVector ydash(halfz * std::tan(ParTheta) * std::cos(ParPhi),
                          halfz * std::tan(ParTheta) * std::sin(ParPhi),
                          halfz / std::cos(ParPhi));
      ydash = ydash.unit();
      G4ThreeVector xdash = Roty.cross(ydash);
      G4ThreeVector zdash = xdash.cross(ydash);
      SideRefVec1 = xdash.unit();
      SideRefVec2 = -ydash.unit();
      SideRefVec3 = -zdash.unit();
    }
    else if (testrand >= Probs[2] && testrand < Probs[3])
    {
      // side is -y
      x = x + z * std::tan(ParTheta) * std::cos(ParPhi) - halfy * std::tan(ParAlpha);
      y = -halfy + z * std::tan(ParTheta) * std::sin(ParPhi);
      // z = z;
      // Cosine-law
      G4ThreeVector ydash(halfz * std::tan(ParTheta) * std::cos(ParPhi),
                          halfz * std::tan(ParTheta) * std::sin(ParPhi),
                          halfz / std::cos(ParPhi));
      ydash = ydash.unit();
      G4ThreeVector xdash = Roty.cross(ydash);
      G4ThreeVector zdash = xdash.cross(ydash);
      SideRefVec1 = xdash.unit();
      SideRefVec2 = ydash.unit();
      SideRefVec3 = zdash.unit();
    }
    else if (testrand >= Probs[3] && testrand < Probs[4])
    {
      // side is +z
      z = halfz;
      y = y + halfz * std::sin(ParPhi) * std::tan(ParTheta);
      x = x + halfz * std::cos(ParPhi) * std::tan(ParTheta) + y * std::tan(ParAlpha);
      // Cosine-law
      SideRefVec1 = Rotx;
      SideRefVec2 = Roty;
      SideRefVec3 = Rotz;
    }
    else if (testrand >= Probs[4] && testrand < Probs[5])
    {
      // side is -z
      z = -halfz;
      y = y - halfz * std::sin(ParPhi) * std::tan(ParTheta);
      x = x - halfz * std::cos(ParPhi) * std::tan(ParTheta) + y * std::tan(ParAlpha);
      // Cosine-law
      SideRefVec1 = Rotx;
      SideRefVec2 = -Roty;
      SideRefVec3 = -Rotz;
    }
    else
    {
      G4cout << "Error: testrand out of range" << G4endl;
      if (verbosityLevel >= 1)
        G4cout << "testrand=" << testrand << " Probs[5]=" << Probs[5] << G4endl;
    }

    RandPos.setX(x);
    RandPos.setY(y);
    RandPos.setZ(z);
  }

  // Apply Rotation Matrix
  // x * Rotx, y * Roty and z * Rotz
  if (verbosityLevel == 2)
    G4cout << "Raw position " << RandPos << G4endl;

  x = (RandPos.x() * Rotx.x()) + (RandPos.y() * Roty.x()) + (RandPos.z() * Rotz.x());
  y = (RandPos.x() * Rotx.y()) + (RandPos.y() * Roty.y()) + (RandPos.z() * Rotz.y());
  z = (RandPos.x() * Rotx.z()) + (RandPos.y() * Roty.z()) + (RandPos.z() * Rotz.z());

  RandPos.setX(x);
  RandPos.setY(y);
  RandPos.setZ(z);

  // Translate
  particle_position = CentreCoords + RandPos;

  if (verbosityLevel >= 1)
  {
    if (verbosityLevel == 2)
      G4cout << "Rotated position " << RandPos << G4endl;
    G4cout << "Rotated and translated position " << particle_position << G4endl;
  }
  if (verbosityLevel == 2)
  {
    G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
  }
}

void ExSPSPosDistribution::GeneratePointsInVolume()
{
  G4ThreeVector RandPos;
  G4double tempx, tempy, tempz;
  G4double x, y, z;
  x = y = z = 0.;
  if (SourcePosType != "Volume" && verbosityLevel >= 1)
    G4cout << "Error SourcePosType not Volume" << G4endl;
  // Private method to create points in a volume
  if (Shape == "Sphere")
  {
    x = Radius * 2.;
    y = Radius * 2.;
    z = Radius * 2.;
    while (((x * x) + (y * y) + (z * z)) > (Radius * Radius))
    {
      x = posRndm->GenRandX();
      y = posRndm->GenRandY();
      z = posRndm->GenRandZ();

      x = (x * 2. * Radius) - Radius;
      y = (y * 2. * Radius) - Radius;
      z = (z * 2. * Radius) - Radius;
    }
  }
  else if (Shape == "SphericalShell")
  {
    G4double RRndm = pow(G4UniformRand() * (pow(Radius, 3) - pow(Radius0, 3)) + pow(Radius0, 3), 1. / 3);
    G4double ThetaRndm = acos(2 * G4UniformRand() - 1);
    G4double PhiRndm = 2 * pi * G4UniformRand();
    x = RRndm * sin(ThetaRndm) * cos(PhiRndm);
    y = RRndm * sin(ThetaRndm) * sin(PhiRndm);
    z = RRndm * cos(ThetaRndm);
  }
  else if (Shape == "Ellipsoid")
  {
    G4double temp;
    temp = 100.;
    while (temp > 1.)
    {
      x = posRndm->GenRandX();
      y = posRndm->GenRandY();
      z = posRndm->GenRandZ();

      x = (x * 2. * halfx) - halfx;
      y = (y * 2. * halfy) - halfy;
      z = (z * 2. * halfz) - halfz;

      temp = ((x * x) / (halfx * halfx)) + ((y * y) / (halfy * halfy)) + ((z * z) / (halfz * halfz));
    }
  }
  else if (Shape == "Cylinder")
  {
    x = Radius * 2.;
    y = Radius * 2.;
    while (((x * x) + (y * y)) > (Radius * Radius))
    {
      x = posRndm->GenRandX();
      y = posRndm->GenRandY();
      z = posRndm->GenRandZ();

      x = (x * 2. * Radius) - Radius;
      y = (y * 2. * Radius) - Radius;
      z = (z * 2. * halfz) - halfz;
    }
  }
  else if (Shape == "Para")
  {
    x = posRndm->GenRandX();
    y = posRndm->GenRandY();
    z = posRndm->GenRandZ();
    x = (x * 2. * halfx) - halfx;
    y = (y * 2. * halfy) - halfy;
    z = (z * 2. * halfz) - halfz;
    x = x + z * std::tan(ParTheta) * std::cos(ParPhi) + y * std::tan(ParAlpha);
    y = y + z * std::tan(ParTheta) * std::sin(ParPhi);
    // z = z;
  }
  else
    G4cout << "Error: Volume Shape Doesnt Exist" << G4endl;

  RandPos.setX(x);
  RandPos.setY(y);
  RandPos.setZ(z);

  // Apply Rotation Matrix
  // x * Rotx, y * Roty and z * Rotz
  tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
  tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
  tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());

  RandPos.setX(tempx);
  RandPos.setY(tempy);
  RandPos.setZ(tempz);

  // Translate
  particle_position = CentreCoords + RandPos;

  if (verbosityLevel == 2)
  {
    G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
    G4cout << "Rotated position " << RandPos << G4endl;
  }
  if (verbosityLevel >= 1)
    G4cout << "Rotated and translated position " << particle_position << G4endl;

  // Cosine-law (not a good idea to use this here)
  G4ThreeVector zdash(tempx, tempy, tempz);
  zdash = zdash.unit();
  G4ThreeVector xdash = Rotz.cross(zdash);
  G4ThreeVector ydash = xdash.cross(zdash);
  SideRefVec1 = xdash.unit();
  SideRefVec2 = ydash.unit();
  SideRefVec3 = zdash.unit();

  if (verbosityLevel == 2)
  {
    G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl;
  }
}

G4bool ExSPSPosDistribution::IsSourceConfined()
{
  // Method to check point is within the volume specified
  if (Confine == false)
    G4cout << "Error: Confine is false" << G4endl;
  G4ThreeVector null(0., 0., 0.);
  G4ThreeVector *ptr;
  ptr = &null;

  // Check particle_position is within VolName, if so true,
  // else false
  G4VPhysicalVolume *theVolume;
  theVolume = gNavigator->LocateGlobalPointAndSetup(particle_position, ptr, true);
  if (!theVolume)
    return (false);
  G4String theVolName = theVolume->GetName();
  G4int theVolID = theVolume->GetCopyNo();
  if (verbosityLevel >= 1)
    G4cout << "volname: " << theVolName << G4endl;
  if (theVolName == VolName)
  {
    if (verbosityLevel >= 1)
      G4cout << "Particle is in volume " << VolName << " CopyID: " << theVolID << G4endl;
    return (true);
    //if(VolID == theVolID || VolID < 0)
    //  return (true);
    //else 
    //  return (false);
  }
  else
    return (false);
}

G4ThreeVector ExSPSPosDistribution::GenerateOne()
{
  //
  G4bool srcconf = false;
  G4int LoopCount = 0;
  while (srcconf == false)
  {
    if (SourcePosType == "Point")
      GeneratePointSource();
    else if (SourcePosType == "Beam")
      GeneratePointsInBeam();
    else if (SourcePosType == "Plane")
      GeneratePointsInPlane();
    else if (SourcePosType == "Surface")
      GeneratePointsOnSurface();
    else if (SourcePosType == "Volume")
      GeneratePointsInVolume();
    else
    {
      G4cout << "Error: SourcePosType undefined" << G4endl;
      G4cout << "Generating point source" << G4endl;
      GeneratePointSource();
    }
    if (Confine == true)
    {
      srcconf = IsSourceConfined();
      // if source in confined srcconf = true terminating the loop
      // if source isnt confined srcconf = false and loop continues
    }
    else if (Confine == false)
      srcconf = true; // terminate loop
    LoopCount++;
    if (LoopCount == 100000)
    {
      G4cout << "*************************************" << G4endl;
      G4cout << "LoopCount = 100000" << G4endl;
      G4cout << "Either the source distribution >> confinement" << G4endl;
      G4cout << "or any confining volume may not overlap with" << G4endl;
      G4cout << "the source distribution or any confining volumes" << G4endl;
      G4cout << "may not exist" << G4endl;
      G4cout << "If you have set confine then this will be ignored" << G4endl;
      G4cout << "for this event." << G4endl;
      G4cout << "*************************************" << G4endl;
      srcconf = true; // Avoids an infinite loop
    }
  }
  return particle_position;
}

void ExSPSPosDistribution::CopyFromSPSP(ExSPSPosDistribution *pos1)
{
  SetPosDisType(pos1->GetPosDisType());
  SetPosDisShape(pos1->GetPosDisShape());
  SetCentreCoords(pos1->GetCentreCoords());
  SetHalfX(pos1->GetHalfX());
  SetHalfY(pos1->GetHalfY());
  SetHalfZ(pos1->GetHalfZ());
  SetRadius(pos1->GetRadius());
  SetRadius0(pos1->GetRadius0());
  ConfineSourceToVolume(pos1->GetConfineVolume());
}
